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A multi-domain spectral method for computing very high precision 3-D stellar models is presented. 
The boundary of each domain is chosen in order to coincide with a physical discontinuity (e.g. 
the star's surface). In addition, a regularization procedure is introduced to deal with the infinite 
derivatives on the boundary that may appear in the density field when stiff equations of state are 
| used. Consequently all the physical fields are smooth functions on each domain and the spectral 

. method is absolutely free of any Gibbs phenomenon, which yields to a very high precision. The power 

y—{ ' of this method is demonstrated by direct comparison with analytical solutions such as MacLaurin 

. „, spheroids and Roche ellipsoids. The relative numerical error reveals to be of the order of 10~ 10 . This 

approach has been developed for the study of relativistic inspiralling binaries. It may be applied to 
a wider class of astrophysical problems such as the study of relativistic rotating stars too. 

! PACS number(s): 02.60.Cb 02.70.Hm 04. 25. Dm 
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One of the most promising source of gravitational waves is the coalescence of inspiralling compact binaries. The 
recent development of interferometric gravitational waves detectors (e.g. GEO600, LIGO, TAMA and VIRGO) gives 
an important motivation for studying this problem. Such a study requires a relativistic formalism to derive the 
equations of motion and then an accurate and tricky method to solve the resulting system of partial differential 
ON ' equations. We have recently Q proposed a relativistic formalism able to tackle the problem of co-rotating as well as 
counter-rotating binaries system, these latter being more relevant from the astrophysical point of view. We present 
O ,■ now a very accurate approach based on multi-domain spectral method that circumvents the Gibbs phenomenon to 
numerically solve this problem and which can be applied to a wide class of other astrophysical situations. 

Various astrophysical applications of spectral methods have been developed in our group (for a review, see Q), 
including 3-D gravitational collapse of stellar core || , neutron star collapse into a black hole , tidal disruption 
of a star near a massive black hole ||. rapidly rotating neutron stars |^-[l2|], magnetized neutron stars [ jl3|Jl4| and 
their resulting gravitational radiation Jl^| , spontaneous symmetry breaking of rapidly rotating neutron stars |16| , |l7| , 
. ,-h | proto- neutron stars evolution p8|-p0[ . 

In computational fluid dynamics, spectral methods are known for their very high accuracy pl| , ^2| : indeed for a C°° 
function, the numerical error decreases as exp(— TV) (evanescent error), where N is the number of coefficients involved 
in the spectral expansion, or equivalently the number of grid points in the physical domain. This is much faster than 
the error decay of finite-difference methods, which behaves as 1/N q , with q generally not larger than 3. For this 
reason, spectral methods are particularly interesting for treating 3-D problems — such as binary configurations — a 
situation in which the number of grid points is still severely limited by the capability of present and next generation 
computers. 

Spectral methods lose much of their accuracy when non-smooth functions are treated because of the so-called Gibbs 
phenomenon. This phenomenon is well known from the most familiar spectral method, namely the theory of Fourier 
series: the Fourier coefficients (c n ) of a function / which is of class C p but not C p+1 decrease as l/n p only. In 
particular, if the function has some discontinuity, its approximation by a Fourier series does not converge towards / 
at the discontinuity point: there remains a gap which is of the order 10%. 

The multi-domain spectral method described in this paper circumvents the Gibbs phenomenon. The basic idea is to 
divide the space into domains chosen so that the physical discontinuities are located onto the boundaries between the 
domains (Sect. |TI|). The simplest example is the case of a perfect fluid star, where two domains may be distinguished: 
the interior and the exterior of the star. The boundary is then simply the surface of the star. The second ingredient of 
the technique is a mapping between the domains defined in this way and some simple mathematical domains, which 
are cross products of intervals: [01,02] x [01,62] x [01,02]. The spectral expansion is then performed with respect to 
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functions of the coordinates spanning these intervals (Sect. III). The method of resolution of a basic equation, namely 
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the Poisson equation, is exposed in Sect. IIIC. For stiff equations of state, the above procedure is not sufficient to 
ensure the smoothness of all the functions. Indeed, for a polytrope with an adiabatic index greater than 2, the density 
field has an infinite derivative on the surface of the star. We present in Sect. IV a method for regularizing the density 
and recover the spectral precision. The power of the multi-domain spectral method is illustrated in Sect. where 
comparisons are performed between numerical solutions obtained by an implementation of the method and analytical 
solutions (ellipsoidal configurations of incompressible fluids). Finally Sect. VI concludes the article by discussing the 
great advantages of the multi-domain spectral method for dealing with relativistic binary neutron stars. 



II. THE PHYSICAL DOMAINS AND THEIR MAPPING 



A. Splitting of the physical space into star-like domains 



In order to treat problems involving Poisson equations with non-compact sources — as they appear in relativistic 
gravitation — we take for the physical domain where the computation must be performed the whole three-dimensional 
space IR 3 . In doing so, we know the physical boundary conditions we have to impose in order to solve the Poisson 
equations. These boundary conditions can be easily set at infinity. We divide IR 3 into Af domains (2?j)o<i<Af-i 
{Af > 2). In the present work, these domains are taken to be star-like (in the mathematical sense) with respect to a 
some origin O, which means that for every point M in the domain T>i, the segment OM is entirely included in {J i<t T>i 
(see Fig. [[]). The multi-domain spectral method we are going to describe can be extended to more general domains, 
at the price of a greater technical (but non conceptual) difficulty. However, for stellar configurations, the star-like 
hypothesis is sufficient for most applications. Let us denote by Si the boundary surface between the domains T>i and 
Pj+i. T>q is simply connected and its boundary is So; we call it the nucleus. For 1 < I < Af — 2, P/'s inner boundary 
is Si-i and outer boundary Si. The last domain, 2?y\/_i, has <Sy\/-2 as inner boundary and extends to infinity (cf. 
Fig. 0). 




FIG. 1. Splitting of the physical three-dimensional space into domains X?o, fi, — ,2?A/"— l (on the figure Af — 3), which are 
star-like with respect to some origin O. The last domain (here D2) extends up to infinity. 

Let us choose some Cartesian frame of IR 3 centered at O and let us call (r, 9, tp) the associated spherical coordinates: 
r G [0, +00 [, G [0,7r] and <p G [0, 2ir[. The mapping of each domain onto the cross product of intervals [01,0,2] x 
[61,62] x [ci,C2] will be defined with respect to the spherical coordinates (r,8,ip). Since each domain T>i is star-like 
with respect to O, the equation of the boundaries Si can be written under the form 

r = Si(6,<p), (1) 

where Si is a smooth function on [0, n] x [0, 27r[. 



B. Mapping of the nucleus 



The basic idea is to introduce a mapping 
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[0,1] x [0,tt] x [0,2tt[ — » 2? 

so that the origin corresponds to £ = and the boundary Sq to £ = const = 1. Using the fact that X>o is star-like, 
a simple form of the mapping (||) can be chosen as 

r = iio(£,0V) (3) 
6» = 0' (4) 

f = <f' , (5) 

where Rq is a smooth function subject to regularity properties which are discussed below. Thanks to Eqs. (0)-(|j|), 
we will make no distinction between 9 and 9', as well as between (p and tp' , i.e. we will abandon the primes on the 
angles. The fact that £>o's boundary coincides with £ = 1 translates as 

RoO-,0,<p)=S o (9,<p) . (6) 

Beside Eq. the function i?o must satisfy some regularity properties induced by the singular behavior of spherical 
coordinates at r — 0, 9 = and 9 = n. We define a function / : 1R 3 — > IR to be regular if it can be expressed as a 
polynomial of the Cartesian coordinates 

x = r sin 9 cos Lp (7) 
y — r sin sin </? (8) 
z = r cos . (9) 

We will assume that all the physical fields are regular functions on each domain T>i (the domains T>i are in fact chosen 
in this manner) with respect to the previous definition. It is easy to see that any regular function is expandable as 

M L 

f(r,e,tp)= ]T sin ' m|0 p *-m ( cos °) eirnv > ( 10 ) 

m—-M e=\m\ 

where L and M are positive integers, L > M, i^_i TO i is some polynomial of degree I — \m\ and T(r 2 ) is some even 
polynomial. 

A simple form of the mapping @-(||) has been already introduced in the literature p3V|25|], namely Ro(£,,8,ip) = 
So(9,tp)!;, where So(9,ip) is the equation of the star's surface [Eq. (Q)]. However for such a mapping the regularity 
condition ( |Io|) would be quite complicated when expressed in terms of (£,8,ip). We choose instead the mapping 
defined by 

Ro(£, 9, <p) = a [f + A (0F (9, <p) + B (£)G (9, <p)] , (11) 

where Aq and Bq are the following even and odd polynomials 

A>(0=3£ 4 -2£ 6 (12) 
B„(0 = (5£ 3 - 3£ 5 )/2 , (13) 

and the constant ckq as well as the two functions Fo(9,ip) and Go(9,ip) are such that (i) the Fourier expansion of 
Fq(9, (p) (resp. Go(9,ip)) with respect to (p contains only even (resp. odd) harmonics and (ii) the equation of the 
surface Sq can be written 

a [1 + F (9, <p) + G Q (9, <p)} = S Q (0, tp) . (14) 
The polynomials Aq and Bq defined by (E3)-(tt3[) are such that 



MO) = B (0) = (15) 

^(0) = B' Q (0) = (16) 

A (l) = B (l) = 1 (17) 

Ml) = B' Q (1) = . (18) 

The properties ( |l4| ) and (|l^) ensure that Eq. (g) is satisfied, i.e. that the mapping ( pi] ) is from [0, 1] x [0, n] x [0, 2tt[ 
to Vq. The polynomials Aq and Bq are chosen in order to satisfy at the minimum level the regularity conditions 
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mentioned above. In particular, near the origin r behaves as r ~ ao£ and is independent of (6,<p), which would not 
have been the case of the mapping r = Sq{9, ip) £ introduced in Refs. |2^-[25[|. 
The Jacobian of the transformation (r, 8, ip) i— > (£. 9, <p) is 

= d(r,9,cp) = dRo 

= a [l + A^)F o (0, ^ + B' o (£)G o (0, <p)] (19) 

Since A' (£) > and B (£) > for any £ € [0, 1], the mapping could become singular if Fq(9, ip) or Gq(9, <p) is negative 
and has an important amplitude. We cannot control the sign of the function Fq(9,p) because it contains only odd 
harmonics of ip. However, by a suitable choice of oiq, one can impose 

G (6,<p)>0, (20) 

as we shall see below. We have found that this condition was sufficient to ensure that J ^ 0, i.e. that the mapping is 
regular, for all the astrophysically relevant situations we have encountered. 

The equation for the surface of the domain X>o being given, under the form (Q), r — Sq(0, tp), the procedure which 
leads to ao, Fq(9, <p) and Gq(8, ip), i.e. to the full determination of the function Rq(£, 8, <fi), is the following one. 

First let us choose a point (0*, ip*) on the surface Sq. Equation (|l4|) implies the following relation 

a [1 + Fo(0„ V*) + Go(0*, ¥>*)] = So(0*,(p.) . (21) 
Let us introduce the following auxiliaries quantities 

lJ,:=ao[Fo(e*,<p*)+Go(0*,<p*)} (22) 
F o (0,tp) :=a F {6,<p) (23) 
Go(0,ip):=aoGo(0,<p)-n. (24) 



Equation ( ^l| ) translates then in 



whereas Eq. (nlh becomes 



a + H = Sq(9*, tp*) , (25) 



F (6,p) + Go(9, v ) = S (9,v) - S o (0*,p*) . (26) 

Having expanded f(6,ip) :— So(9,tp) — Sq(9^,p^,) into Fourier series with respect to ip, one deduces the functions 
F o (0, ip) (resp. G (9, <p)) by taking only the odd (resp. even) harmonics of this Fourier expansion. \i is then computed 
as 



fi = — mm 



{G (fi,<p), (9,tp) e [0,tt] x [0,27r[} . (27) 



In doing so, condition ( |20| ) will automatically be fulfilled. The value of the coefficient ao is deduced from the above 
value of ji via Eq. (Eq). Finally the functions Fq{9, ip) and Gq(9, <p) are computed from Eqs. (|23|) and (p4|). 



C. Mapping of the intermediate domains 

For 1 < I < J\f — 2, we introduce the mapping 



under the form 



[-1,1] x [0,tt] x [0,2tt[ — > 27, ( 

(£,9',<p>) ^ (r,<?,^) W 



r = R l (S,e',<p') (29) 
= 0' (30) 
V = V' , (31) 
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where Ri is a smooth function which satisfies 

R l (-l,9,cp) = Si- 1 (6,ip) (32) 
R l (+l,6,<p) = S l (6,cp) , (33) 

which means that the inner (resp. outer) boundary of T>i is defined by £ = — 1 (rcsp. £ = +1). 
We choose Ri(£,9,(p) as 

R t (£,6,<p) = at [£ + M (£)F (9, <p) + B t (9, ip)} 

+A , (34) 

where Ai and 2?; are the following polynomials 

MO = (f -3£ + 2)/4 (35) 
i? ; (0 = (-? 3 + 3£ + 2)/4, (36) 

and the constants a; and /?/ and the two functions Fi(9, ip) and Gi{9, <p) are defined from the equations of the surfaces 
Si- 1 and Si by 

ai[-l + F l (p,<p)]+p l =Si- 1 (e,tp) (37) 
a, [+1 + Gi{9, ip)] +/3i= Si(9, <p) (38) 
Fi(9,<p)<0 (39) 
Gi(ff, V )>0 (40) 



Note that the polynomials Ai and Bi defined by (|35|)-(36) are such that 

A t {-l) = l and B,(-1) = (41) 

Ai(+1) = and Bi(+1) = 1 (42) 

A[(-l) = Bl(-l) = A[{+1) = = . (43) 

The properties ( p7j ) and ( f4l| ) [resp. d38| ) and p^)] ensure that Eq. (|3^) [resp.^_Eq. (|33|)] is satisfied, i.e. that the 

-i,if x [o,?r] xfo,: 



mapping (|3J) is from [—1,1] x [0, x[0,27r[ to 2?;. The conditi ons (|39| ) and ( |40| ) ensure that this mapping is not 
singular, by the same argument as that presented for Rq in Sect. II B, the sign of Fi(9,<p) being opposite to that of 
Gi(9, <p) because Ai is a decreasing function of £, whereas Bi is an increasing function of £. 

The equation for the inner and outer boundaries of the domain T>i being given, under the form ([!]): r = Si-i(9, if) 
(inner boundary) and r = Si(9,ip) (outer boundary), the procedure which leads to ai, fli, Fi(6,tp) and Gi(9,<p), i.e. 
to the full determination of the function i?z(£, 9, ip), is the following one. 

First let us choose a point (#*,<£>*) on the surface Si~± along with the corresponding point (9*,ip*) on the surface 
Si . Equations (|37| ) and ([38]) imply the following relations 

ai[-l+FK0*,¥>*)]+A=S;-i(0*,¥>*) (44) 
ai[l + G l (p„(p m )]+p l =Si(0.,<p m ) . (45) 

Let us introduce the following auxiliaries quantities 

X:=ouFi(e»,<p») (46) 

fi:=aiGi(e»,tp„) (47) 

F l (9,<p):=a l F l (9,i P )-\ (48) 

Gi(9,i P ):=aiGi(9,i P )- t i. (49) 



Equations ( (44|) and (|45|) then translate in 



-ai + \ + pi = St- 1 (0.,<p m ) (50) 
ai+fi + pi =Si(0*,<p m ) , (51) 



whereas Eqs. ([37]) and (|38|) become 
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Fi(6,<p) = Si_i(ff, ip) - S t -x(0*,<p.) (52) 

G l (e,<p) = Si(e i <p)-Si{o„<p m ) . (53) 

From the values of Fi(9, tp) and Gi(0, <p) obtained above, A and \i are computed as 

A = -max{F i (0,^), (9,<p) G [0,tt] x [0,2tt[} (54) 

/i = -min|G;(^^), (0, <rf G [0, tt] x [0, 2tt[} . (55) 

In doing so, the conditions ( [39| ) and ( |40| ) will automatically be fulfilled. The value of the constants ai and /3; are 
deduced from the above values of A and /i via Eqs. j50| ) and (|5l|). Finally the functions Fi(9,ip) and Gi(6,(p) are 
computed from Eqs. and 



D. Compactification of the infinite domain 

In the case where the external domain T> cxt := X>jv"-i extends to infinity we introduce the mapping 



under the form 



[-1,1] X [0,7T] X [0,27r[ — ► Pcxt 



u:=l/r = U(^e', ( p') (57) 

= 6*' (58) 

= (59) 

where [/ is a smooth function which satisfies 

U-(-l,fl,VJ) = 5™t(fl,VJ)~ 1 (60) 
CT'C+l, 0, v) = , (61) 

where S cx t(9, := Sj\[—2(9, p). The above two equations show that the inner boundary of X> cxt is defined by £ = — 1, 
whereas £ = +1 corresponds to the infinity. We have already introduced such a compactification of the infinite domain 
in Ref. f9J, in the case of a spherical inner boundary. 
We choose the function [/(£, 6, ip) as 

17(6 0, ¥>) = "cxt K + A C xt(OF cx t(0, <p)-l] , (62) 



where A cxt is the same polynomial of £ as that defined in Eq. (J35J), and the constant a ox t and the function F cx ±(9, <p) 
are defined from the equations of the surface iSjvf-2 by 

a cxt [-2 + F«t(0, = SoxtCfl, y)" 1 (63) 
F cxt (0,^)<O. (64) 

The condition (|64| ) ensures that dU/d^ ^ i.e. that the mapping ( |62|) is not singular. 

The equation for the inner boundary of the domain T> ext being given, under the form ([!]): r = Sext(@,<p), the 
procedure which leads to a e xt and F ext (9,ip), i.e. to the full determination of the function U(£,0,<p), is the following 
one. First let us choose a point on the surface <S xt- Equation ( |6^ ) implies the following relation 

a cx t [-2 + F ext (0*, p*)} = ScxtiO*,^*)^ 1 • (65) 
By introducing the auxiliaries quantities 

A := a CKt F cxt (0*, (p*) (66) 
F ext (6, if) := a cxt F cxt (9, <p) - A , (67) 

this equation translates in 
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- 2a ext + A = S eK t(6*,tp*) 1 , (68) 

whereas Eq. ([53]) becomes 

F„t(0,¥>) =^oxt(0,^)" 1 -5' cx t(^,^)" 1 ■ (69) 
From the above value of F cxt (9, ip), A is computed according to 

\ = -min{P ext (9, l p), (9,<p) £ [0,tt] x [0,2tt[} . (70) 

In doing so, the condition (Q) will automatically be fulfilled (recall that a cx t < 0). The value of a e xt is deduced from 
the above value of A via Eq. (|68|). Finally the function F eK t(9, ip) is computed from Eq. j67|). 

III. MULTI-DOMAIN SPECTRAL METHOD 

A. Spectral expansion of a physical field 

The spirit of the multi-domain spectral method is to perform spectral expansions on each domain T>[, and with 
respect to the coordinates (£, 9, ip) instead of the physical coordinates (r,9,<p). We shall take as basis functions 
separable functions of ip), i.e. functions that can be put under the form X(^)0(9)^(ip). The variable ip being 
periodic, it is natural to use Fourier series in tp, i.e. to choose 

$ k (<p) = e lk * -N v /2 <k< N v /2 , (71) 

where AL is an even integer that we will call the number of degrees of freedom in tp. The associated collocation points 
("grid points") are 

tp k = 2tt k/N v < k < N v - 1 . (72) 

As concerns Q{9) one must use functions that are compatible with the expansion ( |To| ) of any regular scalar field /. 
We shall not use sin'" 1 ' 9 Pg-\ m \ (cos 6*), as suggested by Eq. (|l0|), but a wider set, namely the functions 

e kj (9) =cos(j9) 0<j<N g -l form even (73) 
0^(6") =sm(j0) l<j<N e -2 form odd, (74) 

where Ng is an integer that we will call the number of degrees of freedom in 9 and m is the degree of the harmonic 
in the Fourier series with respect to ip: m — k in the present case. The advantages of this choice are to allow the 
use of Fast-Fourier- Transform algorithms for computing the coefficients, as well as very simple matrices for the usual 
differential operators [^6| . The associated collocation points are 

9 J =nj/(N e -l) 0<j<N e -l. (75) 

As concerns the variable £, we also choose a set wider than merely £ , namely: 

Xkji(£)=T2i(Q 0<i<N r -l for j even (76) 
X kji (0=T 2i+1 (0 0<i<N r ~2 for j odd , (77) 

where N r is an integer that we will call the number of degrees of freedom in r and T n denotes the n th degree Chebyshev 
polynomial. The associated collocation points are 

6 = ^(1^) 0<i<N r -l. (78) 

The above choice concerns the nucleus Do only. For the intermediate and external domains, we choose instead 

X kji (0=T i (0, (79) 

along with the collocation points 
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& = - cos(tt i/(N r - 1)) < i < N r - 1 . (80) 

Note that for the nucleus the above choice is the same as that presented in [^7j, once £ is replaced by r. We report 
the interested reader to that paper for a more detailed discussion about this choice (see also Appendices A, B and D 
of II). 

When a symmetry is present, we use different bases, in order to take the symmetry into account. For instance, an 
often existing symmetry is the symmetry with respect to the equatorial plane, i.e. the plane 8 = tt/2. In this case, 
we use, instead of ([73)- (174 



@kj{6) = cos(2j0) form even (81) 
&kj(0) = sin((2j + 1)6») for m odd . (82) 

The associated collocation points span only [0,7r/2], instead of [0,7r]: 

w 3 



'j 



2 N B -1 



< j < N e - 1 . (83) 



Another usual symmetry is the above equatorial symmetry augmented by the symmetry under the transformation 
Lp i — ► ip -\- 7r. This is the case of a triaxial ellipsoid, or of an axisymmetric star perturbed by even m modes. In this 
case, the cp-basis functions are 

$ fc M=e 2 ^. (84) 
The associated collocation points span [0,7r[, instead of [0, 27r[: 

tp k = 7T k/N v < k < N v - 1 . (85) 

The basis in 8 become 

e fci (0)=cos(2j0), (86) 



instead of (81)-(p2[), the collocation points in 8 remaining those given by Eq. (|83|). In this case, the basis for £ in the 
nucleus contain only even polynomials: 

X kji (0 = T«(0 , (87) 
the collocation points remaining those given by Eq. (j78[). 



B. Differential operators 

In this section, we present how a first order differential operator, the gradient, and a second order one, the Laplacian, 
both applied to a scalar field, are expressed in term of the coordinates system described above. The computation of 
any other kind of operator is straightforward. 

The components of the gradient of a scalar field / in an orthonormal basis associated with the spherical coordinates 
(r, 8, ip) are 

df - T- ldf (M) 

^- Jl (88) 

r 88 Ri 88' J x <9£ y ' 

1 df_ 1 df J 3 8f 



r sin 8 dtp Ri sin 8' dp' J\ <9£ 
where the following abbreviations have been introduced: 



J, := ^ (91) 

J - 1 3Rl (02) 

h - RiW (92) 

j ._ 1 dR, 

Ja - R^W ■ ( } 
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Note that we have re-introduced the primes on 9 and ip [cf. Eqs. (§)-(||)] to avoid any confusion between the 
partial derivatives. The partial derivatives that appear in the quantities Ji are computed by (i) a (banded) matrix 
multiplication on the coefficients of the spectral expansion of the functions Fi(0, <p) and Gi(9, (f) and (ii) analytically 
for the polynomials Ai(£) and Bi(£). In the nucleus, J2 is re-expressed as 



■h 



(3£ 3 - 2£ 5 )dF /d9' + ±(5£ 2 - 3^)dG /d9' 
1 + (3£3 - 2?)F Q + (5e - 3£ 4 )Go 



(94) 



in order to avoid any division by a vanishing quantity at £ = 0. The same thing is done for J 3 . 

The above expressions are valid for the nucleus and the intermediate domains, i.e. for I — 0, . . . , J\f — 2. For the 
compactified domain £> C xt, the quantity to be considered is r 2 V/ instead V/. Indeed, gradients in the compactified 
domain are used in the computation of non-linear terms in the relativistic gravitational field equations (scalar products 
of gradients of the metric potentials). We shall see below that the source of the Poisson equation on 2? GX t is to be 
multiplied by r 4 , so that if each gradient is multiplied by r 2 , this multiplication by r 4 is automatically performed. 
The orthonormal components of r 2 V/ on T> ext are 



2 

r x 



r 2 x 



9 

r x 



1 



df 
dr 

15/ 
r d9 

df 



'dU 

U 86' ' 
1 



'9/ 
8U 



df 



r sin 9 dip U sin 9' dtp' 



dU 



1 



dU df 



U sin 9' dtp' d£ 



(95) 
(96) 
(97) 



The Laplacian of a scalar field / reads 



A/ = (1 + fi + H) % + If} + ^1 - *-'{•■ 



d 2 f \ 



-^A 6ip Ri + jf 1 ! Jf x (i + j 2 2 + J 3 2 )^— - -r - 



de 



j 2 dy 

Ri d9d£, ' Ri sin 9 dtpd£ J 
J3 d 2 R 



Rid9d£ Ri sin 9 dtpdS, J 



where the primes on 9 and ip have been abandoned again and the following abbreviation has been introduced: 

d 2 



^ QQ2 + tSin0d0 1 sin 2Qdcp 2 



d 1 

+ 



(98) 



(99) 



C. Resolution of the Poisson equation 



For many astrophysical applications, one has to solve the Poisson-like equation 

A/ = a(f) , 



(100) 



for some 'potential' /. Note that for relativistic computations, cr(f) is not compactly supported (see e.g. and 
generally decreases as 1/r 4 when r — > +00. 

When expressed in terms of the variables j the Laplacian takes the complicated form (|98|), for which it is 

not obvious to find eigcnfunctions. Therefore, we introduce, in each domain T>i, a new radial coordinate 



( := ait; + Pi , 



(101) 

where ai and (3i are the same constants as in Eqs. (|ll]) and ([34]) (in the nucleus: /3q = 0). In the exterior domain, we 
introduce 



V 



«ext(£ - 1) j 



(102) 



where a ex t is the same constant as in Eq. (p2|). We may then split the Laplacian operator A into a pseudo-Laplacian A 
and a part which would vanish if the domains T>i were exactly spherical (in this case, the coordinates ( and r\ introduced 
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here above would coincide with the physical coordinates r and u = 1/r respectively) . By pseudo-Laplacian, we mean 
the operator which once expressed in terms of (C, 0, ip) has the same structure than the Laplacian operator in spherical 
coordinates: 

* 9 2 2 d 1 * 

where Ag v is defined by Eq. (|9^). In the exterior domain, the pseudo-Laplacian is defined instead by 

K ^^f+ ■ (104) 

It is much easier to invert the operator A than the operator A: using spherical harmonics in (0,ip>), the problem 
reduces to a system of second order ordinary differential equations with respect to the variable £. Moreover, the 
junction conditions betw een t he various domains are easily imposed, as explained below. 
The Poisson equation (100) becomes 

aAf = *(J)+H{f), (105) 

where 

a := a 2 J^ 2 (1 + jf + J 2 ) , (106) 



Ri 



(1 + J| + J 2 ) - 1 

J3 



2 5/ 



A" 



;(1 + J| + J|) 



-^&e v Ri + Jf 1 ^(1 + 4 + 4) 



3j oe 



'hd 2 Ri J 3 d 2 R t \ 
.Rid8d£ Ri sin dd(pd£J 



dj_ 



In order to let appear only the operator A in the left-hand-side of Eq. (105), we introduce 



ai 



max a 



and recast Eq. (UOq) into 



A/ 



l r 

ai 



a(f)+K(f) + ( ai -a)Af 



(107) 



(108) 



(109) 



Since / appears on the right-hand-side of this equation, we solve it by iteration. Besides, we introduce some relaxation 
in the computation of the term A/ in the right-hand-side of Eq. (109). More specifically, we solve at each step of the 
iterative scheme the following equation 



A/ 



.7+1 = ~J 



(110) 



where the index J denotes the step at which the quantities are taken and a J is the following source, computed from 
the value of / at the step J: 



^{^(f) + n(f J ) + (a t - a) [Aa' 7 - 1 + (1 - X)a J - 2 } } 



(111) 



In this expression, A is a relaxation parameter (a typical value is A = 1/2) and lZ(f J ) is to be computed according 
to Eq. (107). For the first step (J = 0), f J , d" 7-1 and a J ~ 2 are set to zero or to their value at a previous step in an 
evolutionary scheme. 
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We have exposed the method of resolution of Eq. (110) el sewhe re [2jj9J. Let us simply mention that we first perform 
a transformation from the bases in (Q, if) described in Sect. Ill A (Chebyshev polynomials in cos0 — Fourier expansion 
in ip) to spherical harmonics YJ n (& , p), by means of a matrix multiplication onto the coefficients of the 9 expansion. 
For each value of (£,m), Eq. (|110|) gives then a second order ordinary differential equation with respect to f, the 
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solution of which amounts to invert a banded matrix. Two solutions of the homogeneous equation (A/ = 0) are then 
added in order to connect the solution and its first derivative across the boundaries between the Af domains. More 
precisely, the global boundary condition, generally / — > when r — > +00, is imposed by setti ng the value of / at the 



exterior boundary of the external domain, which is exactly r = +00 as explained in Sect. II D . The matching between 
the various domains amounts then to the resolution of a simple system of 27V — 1 linear equations for the coefficients 
of the homogeneous solutions to be added in each domain. Note that this matching is performed for each value of 
(l,m). 

IV. REGULARIZATION OF THE SOURCE OF POISSON EQUATION 

A. Description of the method 

The analytical properties of the source of the gravitational field at the boundary of the star depend on the equation 
of state (EOS). For a polytrope of adiabatic index 7 (P cx n 7 ), the matter density n behaves as ffVCT-i) where H 
is the specific enthalpy. Consequently, for 7 > 2, the derivative dn/dH has an infinite value for H = 0, i.e. at the 
surface of the star and dn/dr diverges at surface of the star. For values of 7 < 2 only derivatives of higher order 
diverge (actually there exists some value of 7, e.g. 7 = 4/3, for which all derivatives vanish or have a finite value at 
the surface of the star). 

In a stead y sta te configuration H is Taylor expandable at the neighbourhood of the star's surface (this can be easily 



seen on Eq. (123) below). Therefore H vanishes as r — R(9, tp) where R(9, (p) is the equation of the star's surface and 
n behaves as n ~ [r — R(6, t^)] 1 ^ 7-1 ' (this analysis remains valid even for EOS more general than the polytropic one 
if 7 is defined as 7 = d(ln(P) /dln(H)\jj=o)- Consequently n is generaly not a C°° function. This singular behaviour 
implies that the C 2 truncation error of the spectral approximation is no more evanescent and moreover that Gibbs 
phenomenon is present. This fact is especialy awkward when studying the stability of equilibrium configurations 
or looking for bifurcation points because high accuracy is required. In practice 7 cannot be larger than 3 ||l6| , |l7f . 
Note that in the litterature the potential in spherical coordinates is often computed by expanding the source in 
spherical harmonics Yi m (9,ip) and by computing the radial part with a finite difference method. In this case the 
Gibbs phenomenon will appear in the angular part of the solution. The situation is even worse if the radial part of 
the potential is computed with a spectral method. A method to recover spectral accuracy in such cases is as follows. 

We first introduce a known potential $div such that n^w := A<I>div as the same pathological behaviour as n and 
such that n — n^iv is a regular function (at least more regular then n) and numericaly solve 



where $ rC gu := $ - $div 
Consider for instance 



A$ rogu =n- ?idiv (112) 



<& div = F(£,0,^(l-£ 2 ) (a+2) (H3) 



where a = 1/(7 — 1), F is an arbitrary regular function and £ is a new radial variable such that £ = 1 at the surface 
of the star (see Sect. [n]). It is easy to see that A$div has a term vanishing at the surface as (1 — £) Q (i.e. with the 
same pathological behaviour as n). We have indeed 

A$ div = AF£ 2 (1 - £ 2 )(«+ 2 ) - 4(a + 2)f (1 - £ 2 ) (a+1) c\F 

+(a + 2) [-6(1 - £ 2 )( Q+1 ) + 4(a + 1)£ 2 (1 - £ 2 H F (114) 

where A is the Laplacian computed with respect to (£, 9, (p) (cf. Eq. ( |103| )). 

The choice of the factor of (1 — £)( Q + 2 ) is done in order that $ d iv has the required regularity properties at £ = and 
the required behaviour at the boundary of the star. The choice of F(£, 9, ip) is arbitrary. If we choose for F(£, 9, tp) an 



harmonic function, AF = 0, the first term of the r.h.s. of the Eq. (114) vanishes. This is an advantage because this 
term can be quite large and, consequently, give rise to a large error in computing ^regu- We write $div = J2i m a lm^lmi 
where 

*ta :=£'(! -e) {a+2) Yi m (9,<p) (115) 
and where ai m are some numerical coefficients to be determined. We then obtain ndi v (£, 9, </?) = 

Zl m *lmCl(Z)Yr(0,<p), With 
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C/(C) := (a + 2) [-(4/ + 6)(l-^ 2 )( Q+1 ^ i +4(a + l)e ( ' +2) (l-e 2 ) Q ] • (116) 

We have now to determine the values of the coefficients a\ m which give the most regular function n rcgu := n — n<n v . 
The criterion which seems to give the best results is the following one. 

We expand n and as truncated series of spherical harmonics Yi m (9, ip) and Chebyshev polynomial 

I,L,M 

n(t,0,<p)= J2 ^imT t (OYi m (e,ip) (117) 
and each of the functions C/(£) in a Chebyshev series: 

i 

The value of ai m is computed in such a way that the I th coefficient of the truncated series of n rogu vanishes: 

aim = nii m /Cu . (119) 

By means of the above procedure, we eliminate in n the pathological term vanishing as (1 — £)" but we introduce 
another pathological term oc (1 — £ 2 ) a+1 . However the divergence occurs in a higher order derivative of this term so 
that it has a much weaker effect on the accuracy of the result. The method can be improved by taking 

$ div = 6, <p)(l - i 2 ) a+2 [oi + a 2 (l - £ 2 ) + 03(1 - + ■■■+ Ml - t 2 )** 1 ] • ( 12 0) 



instead of (113). The coefficients are chosen in a such a way that the 1 st , 2 , . . . ,K th derivatives of n rcgu vanish 
at £ = 1. Let us call K the regularization degree of the procedure. 

Note that, since $div and d^div vanish at the surface of the star, the boundary condition one has to impose to solve 
— ?i rogu is the same than that for A$ — n. We want to point out that, the above regularization technique can 
be used, mutatis mutandis, also when a finite difference method is used. 



B. Examples 

Consider two polytropic EOS of adiabatic index 7 = 3 and 7 = 10 with a spherically symmetric distribution of the 
enthalpy H = 1 — £ 2 . The corresponding sources density are n^{^) = (1 — ^ 2 ) 1 / 2 and Uio(£) = ( 1 — ^ 2 ) 1 / 9 . Figure || 
shows the mass distributions n and n regu for various values of the regularization degree K [Eq. ( |12C| )]. Note that in 
the case of T" — 10 the procedure improves considerably the behaviour of the source "n-regu even for K — 1. 



y=3 y=10 




FIG. 2. Original and regularized density profiles for 7 = 3 and 7 = 10 polytropes. The regularized profiled are rescaled to 
take the value 1 at the origin. 
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The method can be tested in the case of 7 = 3 by direct comparison with the analytical solution. In this case the 
gravitational field G = d r $ reads 

G = 4 / (1 - u 2 ) 1/2 u 2 du = -Uarcsinr + r(l - r 2 ) 1/2 - 2r(l - r 2 ) 3/2 l (121) 
r z J Q 8H 

Figure |^ shows the relative C 1 error e on G as a function of the number of degrees of freedom N r . The error e follows 
approximately a power law e oc N~P . The dependence of the exponent (3 with respect to the regularization degree K 
is shown in Figure |[ A value as high as of j3 ~ 17 can be achieved with only K = 6. Note that the relation e = TV"' 3 
is only an approximate law. This means that the error tends to become evanescent when the regularization degree 
increases. 




number of coefficients N r 

FIG. 3. Relative C 1 error e on the gravitational field as a function of the number of degrees of freedom A^ r for different 
regularization degrees K. 

n ' ' ' ' ' n 




1 2 3 4 5 6 

Regularization degree K 

FIG. 4. Dependence of the exponent /3 on the regularization degree K. 
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V. ILLUSTRATIVE APPLICATIONS 



A. 3-D stationary configurations 



In this section, we sketch the general structure of a code for computing single star stationary configurations under 
the influence of rotation and a tidal potential. For simplicity we present only the Newtonian case, the relativistic one 
showing no new qualitative feature but simply involving more equations. 

The equilibrium configuration of a cold star rotating rigidly at the angular velocity Q with respect to some inertial 
frame and embedded in a tidal potential <E>tide is governed by the following three equations^] 



AirGp 



(122) 



H + $ Krs 



-(Sir sm9) 2 



$ti 



idc 



const 



(123) 



P(H) 



(124) 



Equation (122) is the Poisson equation linking the gravitational potential $ gra v to the mass density p. Equation (|123J) 
is the first integral which can be derived from the Euler equation governing the (perfect) fluid velocity under the 
stationarity assu mption; this equation relates the specific enthalpy of the fluid H to the internal and external potentials. 
Finally Eq. (124) is the matter equation of state in the zero-temperature approximation. 

The number of domains used for solving this problem is Af = 3 : one domain £> f° r the star (the nucleus, cf. 
Sect. OB), one inter med iate domain T>\ (cf. Sect. II C), the outer boundary of which is spherical and the external 
domain T>2 (cf. Sect. IIP ). In f act, if one would have to treat only the Newtonian case, one domain would be sufficient 
(i.e. the nucleus) because Eq. (122) has a compact support, which is no more true in the relativistic case. 

A solution is specified by the central value of H (or p), H c say, the value of and the expression of $tidc. These 
quantities being given, the iterative method of resolution is as follows. The Af domains are first taken to be ex actly 
spherical. One starts from a very crude densi ty pro file, for instance p — const in T>o. Solving the Poisson equation (122) 
by m eans of the method presented in Sect. Ill C| leads to the gravitational potential <I>grav Inserting i ts va lue into 
Eq. ( |123| ) give a new profile for the specific enthalpy H (the constant on the right-hand-side of Eq. (123) is fully 
determined by the requirement that H = H c at the centre of the star). The surface of the star being defined by 
H = 0, its equation r = Sq(9, ip) [using the notation of Eq. (|])] is found by searching for the equipotential H = in 
the newly determined H field. This defines a new domain T>q. The corresponding mapping Ro(£,0,ip), i.e. the value 
of the constant a n an d the functions Fo(9,(p) and Go(0,<p) [cf. Eq. (|TT|)] is computed according to the procedure 



described in Sect. 



II B 



The new intermediate domain T>\ is defined by the new inner boundary So (the surface of the 
star) and the unchanged spherical o uter boundary Si. The corresponding mapping -Ri(£, 9, ip) is computed according 
to the procedure described in Sect. [IC. The external domain 2?2 remains unchanged. 

The physical location (r,9, ip) of the collocation points Q,£i,9j,ipk) (where I is the domain index) corresponding 
to these new mappings is a priori different from that of the previous mappings, where all the fields were known. 
Therefore, one has to compute the values of the fields at the new collocation points. In the present case, it is sufficient 
to do so only for the specific enthalpy H. In the domain no. I, the collocation point (£$, 9j, ifk) has the physical radial 
coordinate 



r = R((^9j,ip k ) 



(125) 



where the superscript J refers to the step in the iterative procedure: Rf (£, 9, ip) is the current value of the mapping 
of the domain T>i, whereas i? i J ~ 1 (^, 9, ip) is the previous value. Let us denote by (i J_1 (r, 9, ip), E J ~ 1 (r, 9, ip)) the 
inverse mapping at the previous step. This inverse mapping is computed by searching for the zero of the function 
(/,£) i — ► t Ri(£,, 9, ip). The values of H at the collocation points of the new mapping are given by 



H J (l, 6, 9 j ,i Pk ) = B J - X (L J -\r J , B jt <p k ), S J -V, 0j,<pM,<Pk) 



(126) 



1 see [E8| for a discussion of these equations, including the relativistic case. 
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where r J := Rf(£i,9j,<pk)- The value of H on the right-hand-side is to be taken at a point which a priori does not 
coincide with a collocation p oint in £. It is computed by a direct summation, by means the spectral expansion of H. 
Using the notations of Sect. [II A, it writes 



H(l,t;,9,ip)= E 



k=0 



N„ 



E E Hi kji X kji (0 Q k3 {9) 



3=0 



i=0 



(127) 



where the Hikji are the coefficients of H is the domain no. I. Note that from the computational point of view, this 
summation is the most expensive operation of the method: it scales indeed as (N r NgN v ) 2 . It may be possible to 



replace the whole summation (127) by a truncated one or by some interpolation from the value s of H at the collocation 
points, in order to reduce the computational cost. The main advantage of the summation (127) is that it does not 
introduced any additional error in the method: the right-hand-side of Eq. (127) is the value of H at the specified 
point within the spectral accuracy. 

On ce H is computed at the collocation points of the new mapping by means of Eq. (|l2(j ), the equation of state 
(124) is used to find the values of the mass density p at the collocation points. A new iteration may then begin. 

In all the computations we have made, we have found that this procedure converges. For stationary rotating stars in 
general relativity, a rigorous proof of convergence of such iterative method (except for the re-mapping of the physical 
space at each step) has been given by Schaudt & Pfister pj| . 



B. MacLaurin ellipsoids 



The multi-domain spectral method can handle constant density (incompressible matter) rotating bodies without 
any Gibbs phenomenon. With classical spectral methods, the Gibbs phenomenon would have been very severe since 
the density itself, and not some of its derivatives, is discontinuous across the stellar surface for incompressible fluids. 
This gives us the opportunity to quantify the accuracy of the method since exact analytical solutions are known for 
incompressible bodies: the so-called ellipsoidal figures of equilibrium (see e.g. |30f ). Note that an ellipsoid is not a 
particular case fo r the mapping (pi]): all the coefficients of the expansion of Fo(0,ip) and Go(6,ip) onto the bases 
described in Sect. Ill A are non-zero. In this respect, the ellipsoidal figures constitute a strong test of the method. 



MacLaurin ellipsoid 
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FIG. 5. Logarithm of the relative error of the numerical solution with respect to the number of degrees of freedom in 9 for 
a MacLaurin spheroid at the Jacobi-Dedekind bifurcation point (the number of degrees of freedom in r is iV r = 2N$ — 1). Also 
shown is the error in the verification of the virial theorem. 



For single rigidly rotating objects in the Newtonian regime, the more simple ellipsoidal solutions are constituted 
by the family of MacLaurin spheroids, which are a xisymmctric about their rotation axis. We ha ve co mputed them 
by means of the procedure presented in Sect. V A, setting $tidc = and the equation of state (124) to be simply 
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p = const. The axisymmetry allows to employ N v = 1. The code converges towards ellipsoidal configurations and 
we measure the error by comparing the eccentricity e := yl — (r p /r eq ) 2 (where r p and r eq are respectively the polar 
and equatorial radii) of the numerical solution with that of the analytical solution. The result of this comparison is 
presented in Fig. || for a MacLaurin spheroid located on the MacLaurin sequence at the point where the Jacobi and 
Dcdckind sequences branch off: the eccentricity is e = 0.8127, which corresponds to the ratio r p /r cq = 0.5827. Shown 
in Fig. U is the relative error on the eccentricity as a function of the number of coefficients in the 9 expansion. For 
these calculations, the number of coefficients in the £ expansion in each domain is N r = 2Ng — 1. The straight line 
behaviour of the left side of Fig. || shows that the error is evanescent, i.e. that it decreases as exp(— Ng). For Ng > 20, 
the error saturates at the level of 10~ 12 — 10~ n . This is due to the round-off errors in the computation, which is 
performed with a 15-digit accuracy. It is instructive to compare this result with that obtained with a classical spectral 
method, i.e. with fixed spherical domains, as exposed in For instance, the Fig. ^ can be directly compared with 
Fig. 7 of ||: this latter shows an power-law error decay only (of the type A^ 4 ' 5 ), due to the Gibbs phenomenon 
at the star's surface. Moreover, the error saturates at the level of 10 -5 . Note that this result was obtained with a 
polytropic equation of state (adiabatic index 7 = 2), for which the density is continuous across the surface of star; 
the fixed-spherical-domain spectral method presented in Q was not able to treat incompressible fluid. 

Also shown in Fig. ^| is the relative accuracy with which the 3-D virial theorem is satisfied. The 3-D Newtonian 
virial theorem^ states that for a stationary configuration 2T + 3P + W = 0, where T is the total kinetic energy (with 
respect to the inertial frame) , P is the integral of the pressure throughout the star and W is the gravitational potential 
energy. We have computed each of these three integrals for the numerical solution and evaluated the quantity 



2T + 3P 



\W\ 



(128) 



For an exact solution, e = 0. The triangles plotted in Fig. || depict the value of log 10 e for the numerical models. 
Figure [| shows that the virial error is very well correlated with the error evaluated by a direct comparison with the 
analytical solution. This gives us a strong confidence when using the virial error to evaluate the numerical error in 
more general cases, when no analytical solution is available. 

C. Roche ellipsoids 

Roche ellipsoids are equilibrium solutions for incompressible fluid bodies in a synchronized binary system, within 
the approximation of taking only the second order term in the expansion (around the center of mass of one star) of 
the gravitational potential of the companion. They are obtained by setting 

GM comp / x 2x 2 - y 2 - z 2 \ . . 

$tidc = 1 + - + 129 



in Eq. (123), where a is the abscissa of the center of mass of the companion in the Cartesian frame {x, y, z) centered 



at the center of mass of the star under consideration. Note that in Eq. (123), r must now be the distance to the center 
of mass of the binary system and 6 the angle with respect to the rotation axis of the system. Moreover, fi must be 
chosen so that D, 2 = G(M + M comp )/a 3 , in order that the linear term in x which appear in Eq. ( |123| ) vanishes and 
one is left with an ellipsoidal solution. 

The analytical solutions for Roche ellipsoid are given in the classical book by Chandrasekhar |50| . However they are 
given with an accuracy of five digit only (Table XVI in Ref. [^0[), which is not sufficient for our comparison project: 
the accuracy achieved by the numerical code is far better than 10 -5 as we shall see below. Therefore, we have written 
a small Mathematica [^|| program to compute Chandrasekhar's "index symbol" A\, A2 and A3 and obtain Roche 
solutions with an arbitrary number of digits. 

Figure ^ present the results of the comparison between the numerical solution obtained by means of the method 



described in Sect. V A and the analytical solution. Let us recall that ellipsoidal shapes are not privileged in our 
formalism, so that this type of comparison constitute a strong test of our method. The comparison is conducted at 
fixed values of fl/p and the mass ratio M comp /M. Two global errors can be then defined: (i) the error on the axis ratio 
02/01 and (ii) the error on the axis ratio 03/01, a\ being the major axis of the triaxial ellipsoid (directed along the 
line of the two centers of mass), 02 being the orthogonal axis in the orbital plane, and 03 being the axis perpendicular 



2 as opposed to the 2-D virial identity, see |Slj and |32| for a discussion 
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to the orbital plane. These two errors are shown in Figure ^| for a Roche ellipsoid with fi 2 /{irGp) = 0.1147 and 
M comp /M = 1. The corresponding axis ratios are aija\ = 0.7506 and 03/01 = 0.6853. The numerical solution is 
depicted in Fig. 0by three plane sections obtained with the following (small) numbers of coefficients: N r — 13, Ng = 7 
and N v — 6. Also show in this Figure is the numerical grid (collocation points) used in the problem (only the domain 
2?o and a part of X>i are represented in the Figure). Even with such a small number of points, the relative error is of 
order 1 x 10~ 4 (cf. Fig. ||) ! This explains why despite the numerical grid is quite coarse, the iso-enthalpy surfaces 
shown in Fig. |^ are so smooth. 

Figure ^| gives the two global errors as a functions of the number of coefficients in the 9 expansions, Ng . The number 
of coefficients employed in the other directions are N r = 2Ng — 1 and N v = Ng — 1. As in Fig. [| the exponential 
decay of the error for Ng < 13 means that the error is evanescent. For Ng > 19, the error saturates at the level of a 
few 10 -10 due to the round-off errors in the computation, this latter being performed with a 15-digit accuracy. The 
cost in CPU time for different numbers of degrees of freedom is shown in Table |. 

Roche ellipsoid 



a 2 /a, = 0.75 a 3 /a, = 0.68 
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FIG. 6. Logarithm of the relative global error of the numerical solution with respect to the number of degrees of freedom 
in for a Roche ellipsoid for an equal mass binary system and fi 2 /(irGp) = 0.1147 (the numbers of degrees of freedom in the 
other directions are N r = 2Ng — 1 and N v = N$ — 1). Also shown is the error in the verification of the virial theorem. 



TABLE I. CPU time cost on a R4400/150MHz processor as a function of the number of degrees of freedom for the calculation 
of the Roche ellipsoid configuration corresponding to Fig. ^j. The iteration is halted when the relative discrepancy between two 
successive steps reaches 10~ 13 . 
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FIG. 7. Orthogonal plane sections in the numerical solution obtained for the Roche ellipsoid represented by the second set 
of points starting from the left on Fig. ^ (i.e. corresponding to N r — 13, N$ = 7 and N v = 6). Shown are the iso-enthalpy 
lines, as well as the numerical grid. This computation took a few seconds on a R4400/150MHz processor. 



Also shown in Fig. |g is the relat ive a ccuracy with which the 3-D virial theorem is satisfied. This error estimator is 
defined in the same way as in Sect. VB. As in the axisymmetric case (MacLaurin ellipsoids), we find a high correlation 
between the virial error and the errors obtained by direct comparison with the analytical solution. 



VI. CONCLUSION AND PERSPECTIVES 



We have presented a new numerical approach capable of handling the surface discontinuities of stellar configurations, 
provided these discontinuities are star-like, which covers a wide range of astrophysically relevant situations. When used 
along with spectral methods this adaptive-domain technique ensures that no Gibbs phenomenon can appear. This 
results in a very high precision (evanescent error), as demonstrated in Sect. |v| by comparison with exact analytical 
solutions. The relative error for 3-D configurations can reach 10~ 10 with a relatively small number of degrees of 
freedom (N r x No x N v = 37xl9xl8in each domain). Let us recall that very high accuracy is required for a lot 
of astrophysical problems such as numerical stability analysis. Among these problems let us mention the study of 
symmetry breaking of rapidly rotating stars and the determination of the orbital frequency of the last stable orbit of 
a neutron stars binary system. 

The multi-domain spectral method is particularly well adapted to the computation of relativistic binary neutron star 
system. Three sets of domains can be used in this problem (see Fig. ||): two sets of (three or more) domains centered 
on each star and a third set of (two or more) domains centered at the intersection between the rotation axis and the 
orbital plane. This latter set of domains which reaches spatial infinity is required to compute the gravitational field 
of relativistic configurations. When needed, the quantities computed on one of the three domain sets are evaluated at 
the collocation points of another set by means of the method presented in Sect. V A. We are currently applying this 



numerical method to the computation of steady-state configurations of relativistic counter-rotating (i.e. irrotational 
with respect to an inertial frame) neutron star binaries, following the formulation developed in |IJ . We will report on 
the astrophysical results in a forthcoming paper. 

An interesting byproduct of the present technical paper is the following one. In a previous work B, we had been able 
to demonstrate that the virial error is representative of the true error (measured by direct comparison with analytical 
solutions) only in the spherically symmetric case. We had inferred that this remains valid in the axisymmetric and 
3-D cases. In the present work, we have confirmed this conjecture, thanks to the ability of the present method to 
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treat incompressible fluids, for which 3-D analytical solutions are available. 




FIG. 8. Representation of the numerical domains that we use to compute relativistic steady-state configurations of binary 
neutron stars systems. The external domain extends to spatial infinity in order to compute the exact gravitational potentials. 
Due to the symmetry of the problem, only the z > part of space is taken into account. 
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